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Abstract 

A brief introduction to the projector augmented wave method is given 
and recent developments are reviewed. The projector augmented wave 
method is an all-electron method for efficient ab-initio molecular dynam- 
ics simulations with full wave functions. It extends and combines the 
traditions of existing augmented wave methods and the pseudopotential 
approach. Without sacrificing efficiency, the PAW method avoids transfer- 
ability problems of the pseudopotential approach and it has been valuable 
to predict properties that depend on the full wave functions. 

1 Introduction 

The main goal of electronic structure methods is to solve the Schrodinger equa- 
tion for the electrons in a molecule or solid, to evaluate the resulting total 
energies, forces, response functions and other quantities of interest. In this pa- 



per we review the Projector augmented wave (PAW) method |B16chl 1994|, an 
electronic structure method for ab-initio molecular dynamics with full wave func- 
tions. The main goal of this paper is not to provide a particularly complete or 
detailed account of the methodology, but rather to lay out the underlying ideas. 
A more rigorous description can be found in the origi nal paper p31ochl 1994 1. 



Density functional theory iHohcnberg et al 1964, Kohn et al 1965 maps a 



description for interacting electrons onto one of non-interacting electrons in an 
effective potential. The remaining one-electron Schrodinger equation still poses 
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substantial numerical difficulties: (1) in the atomic region near the nucleus, the 
kinetic energy of the electrons is large, resulting in rapid oscillations of the wave 
function that require fine grids for an accurate numerical representation. On 
the other hand, the large kinetic energy makes the Schrodingcr equation stiff, 
so that a change of the chemical environment has little effect on the shape of 
the wave function. Therefore, the wave function in the atomic region can be 
represented well already by a small basis set. (2) In the bonding region between 
the atoms the situation is opposite. The kinetic energy is small and the wave 
function is smooth. However, the wave function is flexible and responds strongly 
to the environment. This requires large and nearly complete basis sets. 

Combining these different requirements is non-trivial and various strategies 
have been developed. 

• Most appealing to quantum chemists has been the atomic point of view. 

Basis functions are chosen that resemble atomic orbitals. They exploit 
that the wave function in the atomic region can be described by a few 
basis functions, while the bonding is described by the overlapping tails of 
these atomic orbitals. Most techniques in this class are a compromise of 
a well adapted basis set with complex matrix elements on the one hand 
and on the other hand numerically convenient basis functions such as 
Gaussians, where the inadequacies are compensated by larger basis sets. 

• Pseudopotentials regard an atom as a perturbation of the free electron 
gas. The most natural basis functions are plane waves. Plane waves are 
complete and well adapted to sufficiently smooth wave functions. The 
disadvantage of the large basis sets required is offset by the extreme sim- 
plicity to evaluate matrix elements. Finite plane wave expansions are, 
however, absolutely inadequate to describe the strong oscillations of the 
wave functions near the nucleus. In the pseudopotential approach the 
Pauli repulsion of the core electrons is therefore described by an effective 
potential that expels the valence electrons from the core region. The re- 
sulting wave functions are smooth and can be represented well by plane 
waves. The price to pay is that all information on the charge density and 
wave functions near the nucleus is lost. 

• Augmented wave methods compose their basis functions out of atom-like 

partial waves in the atomic regions and a set of functions, called enve- 
lope functions, appropriate for the bonding in between. Space is divided 
accordingly into atom-centered spheres, defining the atomic region, and 
an interstitial region for the bonds. The partial solutions of the differ- 
ent regions, are matched at the interface between atomic and interstitial 
regions. 

The projector augmented wave method is an extension of augmented wave meth- 
ods and the pseudopotential approach, which combines their traditions into a 
unified electronic structure method. 

After describing the underlying ideas of the various methods let us briefly re- 
view the history of augmented wave methods and the pseudopotential approach. 
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We do not discuss the atomic-orbital based methods, because our focus is the 
PAW method and its ancestors. 

The augmented wave methods have been introduced in 1937 by Slater [ Slater 1937[ 
and later modified by Korringa |Korringa 1947 , Kohn and Rostokker |Kohn et al 1954 1 . 
They approached the electronic structure as a scattered electron problem. Con- 
sider an electron beam, represented by a plane wave, traveling through a solid. 
It undergoes multiple scattering at the atoms. If for some energy, the outgo- 
ing scattered waves interfere destructively, a bound state has been determined. 
This approach can be translated into a basis set method with energy dependent 
and potential dependent basis functions. In order to make the scattered wave 
problem tractable, a model potential had to be chosen: The so-called muffin-tin 
potential approximates the potential by a constant in the interstitial region and 
by a spherically symmetric potential in the atomic region. 

The pseudopotential approach traces back to 1940 when C. Herring in- 
vented the orthogo nalized plane w ave method [Herring 1940 . Later, Phillips 



[Phillips et al 1959 and Antoncik [Antoncik 195£[ replaced the orthogonality 
condition by an effective potential, that compensates the electrostatic attrac- 
tion by the nucleus. In practice, the potential was modified for example by 
cutting off the singular potential of the nucleus at a certain value. This was 
done with a few parameters that have been adjusted to reproduce the measured 
electronic band structure of the corresponding solid. 

Augmented wave and pseudopotential methods reached adulthood in the 
1970s: At first O.K. Andersen jAndersen 1975[ showed that the energy depen- 
dent basis set of Slater's APW method can be mapped onto one with energy 
independent basis functions, by linearizing the partial waves for the atomic re- 
gions in energy. In the original APW approach the zeros of an energy dependent 
matrix had to be determined, which is problematic, if many states lie in a small 
energy region as for complex systems. With the new energy independent ba- 
sis functions, however, the problem is reduced to the much simpler generalized 
eigenvalue problem, which can be solved using efficient numerical techniques. 
Furthermore, the introduction of well defined basis sets paved the way for full- 
potential calculations. In that case the muffin-tin approximation is used solely 
to define the basis set. The matrix elements of the Hamiltonian are evaluated 
with the full potential. 

Hamann, Schliiter and Chiang [ Hamann et al 1979 showed in 1979 how 
pseudopotentials can be constructed in such a way, that their scattering proper- 
ties are identical to that of an atom to first order in energy. These first-principles 
pseudopotentials relieved the calculations from the restrictions of empirical pa- 
rameters. Highly accurate calculations have become possible. A main disad- 
vantage of these pseudopotentials has been the large basis set size required 
especially for first-row and transition metal atoms. 

In 1985 R. Car and M. Parrinello published the ab-initio molecular dynamics 
method [Car et al 1985 [. Simulations of the atomic motion have become possi- 
ble on the basis of state-of-the-art electronic structure methods. Besides mak- 
ing dynamical phenomena and finite temperature effects accessible to electronic 
structure calculations, the ab-initio molecular dynamics method also introduced 
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a radically new way of thinking into electronic structure methods. Diagonaliza- 
tion of a Hamilton matrix has been replaced by classical equations of motion for 
the wave function coefficients. If one applies friction, the system is quenched to 
the ground state. Without friction truly dynamical simulations of the atomic 
structure are performed. Electronic wave functions and atomic positions are 
treated on equal footing. 

The Car-Parrinello method had been implemented first for the pseudopoten- 
tial approach. There seemed to be unsurmountable barriers against combining 
the new technique with augmented wave methods. The main problem was re- 
lated to the potential dependent basis set used sofar: the Car-Parrinello method 
requires a well defined and unique total energy functional of atomic positions 
and basis set coefficients. Therefore, it was one of the main goals of the PAW 
method to introduce energy and potential independent basis sets that were as 
accurate and numerically efficient as the previously used augmented basis sets. 
Other requirements have been: (1) The method should match the efficiency of 
the pseudopotential approach for Car-Parrinello simulations. (2) It should be- 
come an exact theory when converged and (3) its convergence should be easily 
controlled. We believe that these criteria have been met, which explains why 
the PAW method becomes increasingly wide spread today. 

We would like to point out that most of these seemingly singular devel- 
opments did not come out of the blue, but the ideas seemed to have evolved 
in the community. In the case of the PAW method, similar ideas have been 
developed by Vanderbilt [Vanderbilt 199C| in the context of ultra-soft pseu- 
dopotentials. The first dynamical simulations using a semiempirical electronic 
structure method have been performed by Wang and Karplus | Wang et al 1973| 
in 1973. The first ab-initio pseudopotentials have been published by Zunger 
[Zunger et al 197^ one year before Hamann, Bachelet and Schliiter [Hamann et al 1979 1. 



2 Transformation theory 

At the root of the PAW method lies a transformation, that maps the true wave 
functions with their complete nodal structure onto auxiliary wave functions, that 
are numerically convenient. We aim for smooth auxiliary wave functions, which 
have a rapidly convergent plane wave expansion. With such a transformation 
we can expand the auxiliary wave functions into a convenient basis set, and 
evaluate all physical properties after reconstructing the related physical (true) 
wave functions. 

Let us denote the physical one-particle wave functions as |^'n) and the aux- 
iliary wave functions as jvPn). Note that the tilde refers to the representation 
of smooth auxiliary wave functions, n is the label for a one-particle state and 
contains a band index, a /c-point and a spin index. The transformation from 
the auxiliary to the physical wave functions is T. 

I*n) = T\^n) (1) 

We use here Dirac's Bra and Ket notation. A wave function ^„(r) corre- 
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sponds to a kct the complex conjugate wave function ^'^(r) corresponds 

to a bra {'^n], and a scalar product / d^r^** (r)\l/m(r) is written as (\l/„|^'f„)- 
Vectors in the 3-d coordinate space are indicated by boldfaced symbols. 

The electronic ground state is determined by minimizing a total energy func- 
tional -E[\l/n] of the density functional theory. The one-particle wave functions 
have to be orthogonal. This constraint is implemented with the method of 
Lagrange multipliers. We obtain the ground state wave functions from the ex- 
tremum condition for 

F{[^n], ^m,n) = - ^[(*„|*™) - 6n,m]An,m (2) 

with respect to the wave functions and the Lagrange multipliers An^m- The 
extremum condition for the wave functions has the form 



where the /„ are the occupation numbers and H = — jlr^^ ''^eff(r) is the 
effective onc-particle Hamilton operator. 

After a unitary transformation that diagonalizes the matrix of Lagrange 
multipliers Am,n, we obtain the Kohn-Sham equations. 

iJ|*„) = |vl/„)e„ (4) 

The one-particle energies e„ are the eigenvalues of An,m 2/"^/" • 

Now we express the functional F in terms of our auxiliary wave functions 

F{[T^n],Am,n) = E[T^r,] - ^[(*„|rtr|*„) - <5„,„]A„,„ (5) 

n,m 

The variational principle with respect to the auxiliary wave functions yields 

Ttifr|*„) = rtT|*„)e„. (6) 

Again we obtain a Schrodinger-like equation, but now the Hamilton operator 
has a different form, T^HT, an overlap operator T^T occurs, and the resulting 
auxiliary wave functions are smooth. 

When wc evaluate physical quantities we need to evaluate expectation values 
of an operator A, which can be expressed in terms of either the true or the 
auxiliary wave functions. 

(A) = ^/„(*„|A|*„) = ^/„(*„|rtAr|*„) (7) 

n n 

In the representation of auxiliary wave functions we need to use transformed 

operators A = T^AT. As it is, this equation only holds for the valence electrons. 
The core electrons are treated differently as will be shown below. 
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The transformation takes us conceptionally from the world of pseudopo- 
tentials to that of augmented wave methods, which deal with the full wave 
functions. We will see that our auxiliary wave functions, which are simply the 
plane wave parts of the full wave functions, translate into the wave functions of 
the pseudopotential approach. In the PAW method the auxiliary wave functions 
are used to construct the true wave functions and the total energy functional 
is evaluated from the latter. Thus it provides the missing link between aug- 
mented wave methods and the pseudopotential method, which can be derived 
as a well-defined approxi mation of the PAW method. 

In the original paper [Blochl 1994 1, the auxiliary wave functions have been 
termed pseudo wave functions and the true wave functions have been termed 
all-electron wave functions, in order to make the connection more evident. I 
avoid this notation here, because it resulted in confusion in cases, where the 
correspondence is not clear cut. 



3 Transformation operator 

Sofar, we have described how we can determine the auxiliary wave functions of 
the ground state and how to obtain physical information from them. What is 
missing, is a definition of the transformation operator T. 

The operator T has to modify the smooth auxiliary wave function in each 
atomic region, so that the resulting wave function has the correct nodal struc- 
ture. Therefore, it makes sense to write the transformation as identity plus a 
sum of atomic contributions Sr 

T=l + Y,Sn. (8) 

R 

For every atom, Sr adds the difference between the true and the auxiliary wave 
function. The index i? is a label for an atomic site. 

The local terms Sr are defined in terms of solutions of the Schrodinger 
equation for the isolated atoms. This set of partial waves will serve as a 
basis set so that, near the nucleus, all relevant valence wave functions can be 
expressed as superposition of the partial waves with yet unknown coefficients. 

*(r) = J2 Mr)cr for |r - Rr\ < rc,R (9) 

ieR 

The index i refers to a site index _R, the angular momentum indices {£, m) and an 
additional index that differentiates partial waves with same angular momentum 
quantum numbers on the same site. With i G R we indicate those partial waves 
that belong to site R. Rr is the position of the nucleus of site R. 

Note that the partial waves are not necessarily bound states and are therefore 
not normalizable, unless we truncate them beyond a certain radius r^fl. The 
PAW method is formulated such that the final results do not depend on the 
location where the partial waves are truncated, as long as this is not done too 
close to the nucleus. 
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Since the core wave functions do not spread out into the neighboring atoms, 
we will treat them differently. Currently we use the frozen-core approximation 
so that density and energy of the core electrons are identical to those of the 
corresponding isolated atoms. The transformation T shall produce only wave 
functions orthogonal to the core electrons, while the core electrons are treated 
separately. Therefore, the set of atomic partial waves includes only valence 
states that are orthogonal to the core wave functions of the atom. 

For each of the partial waves we choose an auxiliary partial wave \(t>i). The 
identity 

\^^) = {1 + SrM,) for leR 
Sum - 10.) -10.) (10) 

defines the local contribution Sr to the transformation operator. Since 1 + Sr 
shall change the wave function only locally, we require that the partial waves 
and their auxiliary counter parts \(f)i) are pairwise identical beyond a certain 
radius re- 

(l)i{r) ^ cf>i{r) for ieR and |r - R/?] > rc^ij (11) 

In order to be able to apply the transformation operator to an arbitrary aux- 
iliary wave function, we need to be able to expand the auxiliary wave function 
locally into the auxiliary partial waves. 

^(r) = ^0,(r)(p,|^) for |r - R^l < r,,^ (12) 

which defines the projector functions The projector functions probe the 

local character of the auxiliary wave function in the atomic region. Examples of 
projector functions are shown in Fig. ^ From Eq. |l^ we can derive J2i \4>i){Pi \ = 
1, which is valid within r^. It can be shown by insertion, that the identity 
Eq. |l^ holds for any auxiliary wave function |^) that can be expanded locally 
into auxiliary partial waves \4>i), if 

{p,\4>j) ^ S,,j for iJeR (13) 

Note that neither the projector functions nor the partial waves need to be or- 
thogonal among themselves. 

By combining Eq. ^ and Eq. 12, we can apply Sr to any auxiliary wave 
function. 

Sr\^) = ^5«|0,)(ft|^')=^(|4)-|0,))(K|§) (14) 

i£R ieR 

Hence the transformation operator is 

T=l + J2[\cl>,)~\4>^,)){p,\ (15) 
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Figure 1: Top: projector functions of the CI atom for two s-type partial waves, 
middle: p-type, bottom: d-type. 



where the sum runs over all partial waves of all atoms. The true wave function 
can be expressed as 

I*) = 1^) +E(l'^«) - l^^))(P.I*> = l^n) - (16) 

i R 

with 

\^r) = Y.\^^){P^m (17) 

ieR 

\^r) = T.\^^){P^\^) (18) 
i£R 

In Fig. H the decomposition of Eq. |l^ is shown for the example of the bonding 
p-cr state of the CI2 molecule. 

To understand the expression for the true wave function, Eq. let us con- 
centrate on different regions in space. (1) Far from the atoms, the partial waves 
are, according to Eq. |ll|, pairwise identical so that the auxiliary wave function is 
identical to the true wave function ^'(r) = ^(r). (2) Close to an atom, however, 
the true wave function ^'(r) = 5'}j.(r) is built up from partial waves that contain 
the proper nodal structure, because the auxiliary wave function and its partial 



wave expansion are equal according to Eq. 12 



In practice the partial wave expansions are truncated. Therefore, the identity 



of Eq. 12 does not hold strictly. As a result the plane waves also contribute to 
the true wave function inside the atomic region. This has the advantage that 
the missing terms in a truncated partial wave expansion are partly accounted 
for by plane waves, which explains the rapid convergence of the partial wave 
expansions. 

Frequently, the question comes up, whether the transformation Eq. |l^ of the 
auxiliary wave functions indeed provides the true wave function. The transfor- 
mation should be considered merely as a change of representation analogous to a 
coordinate transform. If the total energy functional is transformed consistently, 
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Figure 2: Bonding p-cr orbital of the CI2 molecule and its decomposition of the 
wave function into auxiliary wave function and the two one-center expansions. 
Top-left: True and auxiliary wave function; top-right: auxiliary wave function 
and its partial wave expansion; bottom-left: the two partial wave expansions; 
bottom-right: true wave function and its partial wave expansion. 

its minimum will yield an auxiliary wave function that produces a correct wave 
function j^*). 

4 Expectation values 

Expectation values can be obtained either from the reconstructed true wave 
functions or directly from the auxiliary wave functions 

(A) = ^/„(VI/„|A|*„) + E('^"l^l'^n) 
n n—1 

AT, 

= ^/„(*„|rt^r|§„)-f ^(0^1^10$=,) (19) 

n n—1 

where /„ are the occupations of the valence states and Nc is the number of core 
states. The first sum runs over the valence states, and second over the core 
states \4>n)- 

Now we can decompose the matrix elements into their individual contribu- 
tions according to Eq. 

R B.' 
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part 1 



part 2 

Y,{'^R~'^R\A\^R'~'fR') (20) 



part 3 



Only the first part of Eq. 20 , is evaluated explicitly, while the second and third 
parts of Eq. |2^ are neglected, because they vanish for sufhciently local opera- 
tors as long as the partial wave expansion is converged: The function ^E'Jj — 
vanishes per construction beyond some augmentation region, because the par- 
tial waves are pairwise identical beyond that region. The function 5" — 
vanishes inside the augmentation region, if the partial wave expansion is suffi- 
ciently converged. In no region of space both fmictions '^'r — '^r and ^ — are 
simultaneously nonzero. Similarly the functions '^]^ — from different sites 
are never non-zero in the same region in space. Hence, the second and third 
parts of Eq. |2^ vanish for operators such as the kinetic energy ^^V^ and the 
real space projection operator |7')(r|, which produces the electron density. For 
truly nonlocal operators the second and third parts of Eq. ^ would have to be 
considered explicitly. 

The expression, Eq. for the expectation value can therefore be written 

as 

(A) = 5]/„((*„|A|*„) + {K\A\K) + E('^"l^l'^«) 

n n=l 

= J2fn{^n\A\^n)+J2{k\A\k) 
n 71—1 

+ E( E + E i^nim.)) 

R. i.jeR neR 

- E( E D.A4>,\m) + E (21) 

R iyjeR n£R 

where Dij is the one-center density matrix defined as 

A,, = E U^n\Pj){pA^n) = E(P'I*")/"(*"IW) (22) 
n n 

The auxiliary core states, allow to incorporate the tails of the core wave 
function into the plane wave part, and therefore assure, that the integrations of 
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partial wave contributions cancel strictly beyond re- They are identical to the 
true core states in the tails, but are a smooth continuation inside the atomic 
sphere. It is not required that the auxiliary wave functions are normalized. 
For example, the electron density is given by 

n(r) = n(r) + 5](n^(r)-4(r)) (23) 

Ft, 

Mr) = ^/„*:(r)*„(r)+ne 

n 

nR{r) = ^ D,j(j)*j{T)4)i{v) + n^^R 

^kW = ^ Aj^*(r)^i(r)+ne,fl (24) 

ijeR. 

where nc,R is the core density of the corresponding atom and nc,R is the auxiliary 
core density that is identical to nc,R outside the atomic region and a smooth 
continuation inside. 

Before we continue, let us discuss a special point: The matrix element of 
a general operator with the auxiliary wave functions may be slowly converging 
with the plane wave expansion, because the operator A may not be well behaved. 
An example for such an operator is the singular electrostatic potential of a 
nucleus. This problem can be alleviated by adding an intelligent zero: If an 
operator B is purely localized within an atomic region, we can use the identity 
between the auxiliary wave function and its own partial wave expansion 

= (*„|B|*„) - (25) 

Now we choose an operator B so that it cancels the problematic behavior of the 
operator A, but is localized in a single atomic region. By adding B to the plane 
wave part and the matrix elements with its one-center expansions, the plane 
wave convergence can be improved without affecting the converged result. 

5 Total Energy 

Like wave functions and expectation values also the total energy can be divided 
into three parts. 

E{[^n],Ri) = E + Y,{E}i-E},) (26) 

R 

The plane-wave part E involves only smooth functions and is evaluated on equi- 
spaced grids in real and reciprocal space. This part is computationally most 
demanding, and is similar to the expressions in the pseudopotential approach. 
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^3 /^3^.[^^(r) + ^(r)][n(r') + ^(r')] 



+ J d3rn(r)e,e(r, [n]) + J d''rv{v)h{v), (27) 

where Z{r) is an angular dependent core-like density that will be described in 
detail below. The remaining parts can be evaluated on radial grids in a spherical 
harmonics expansion. The nodal structure of the wave functions can be properly 
described on a logarithmic radial grid that becomes very fine near nucleus, 



+ " ,^3^. K(r)+Z(r)][n^(rO + Z(rO] 



=2 



+ ydW(r)e,e(r,K]) (28) 



Sttco j j |r — r'l 

+ j d^rh\r)€^e{r,[n^]) + J d^rv{r)fi\r) (29) 

The nuclear charge density —eZ(r) is defined as a sum of (5-functions on the 
nuclear sites, Z(r) — — '^jiZ]i5{r — R), with the atomic numbers Zj^. Note 
that the self energy of a point charge is infinite and must be subtracted out. 

The compensation density Z{r) — ■2^_R(r) is given as a sum of angular 
momentum dependent Gauss functions, which have an analytical Fourier trans- 
form. A similar term occurs also in the pseudopotential approach. In contrast 
to the norm-conserving pseudopotential approach however, the compensation 
charge is non-spherical and it is constantly adapting to the instantaneous envi- 
ronment. It is constructed such that the augmentation charge densities 

ni^(r)+Z«(r)-fi^(r)-Z^(r) (30) 

have vanishing electrostatic multi-pole moments for each atomic site. As a 
result the sum of all one-center contributions from one atom does not produce 
an electrostatic potential outside their own atomic region. This is the reason 
that the electrostatic interaction of the one-center parts between different sites 
vanish. 

The compensation charge density as given here is still localized within the 
atomic regions, but a technique similar to an Ewald summation allows to replace 
it by a very extended charge density. Thus we can achieve, that all functions in 
E converge as fast as the auxiliary density itself. 
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The potential v, which occurs in Eqs. 27 and 
the form of a zero described in Eq. 



enters the total energy in 

n i,j 

The main reason for introducing this potential is that the self-consistent poten- 
tial resulting from the plane wave part is not necessarily optimally smooth. The 
potential v allows to influence the plane wave convergence beneficially, without 
changing the converged result, v must be localized within the augmentation 
region, where equation O holds. 



6 Approximations 

Once the total energy functional provided in the previous section has been 
defined, everything else follows: Forces are partial derivatives with respect to 
atomic positions. The potential is the derivative of the potential energy with re- 
spect to the density, and the Hamiltonian follows from derivatives with respect 
to wave functio ns. The fictitious Lagrangian approach of Car and Parrinello 



[ Car et al 1985 does not allow any freedom in the way these derivatives are 
obtained. Anything else than analytic derivatives will violate energy conserva- 
tion in a dynamical simulation. Since the expressions are straightforward, even 
though rather involved, wc will not discuss them here. 

All approximations are incorporated already in the total energy functional 
of the PAW method. What are those approximations? 

• Firstly we use the frozen core approximation. In principle this approxi- 
mation can be overcome. 

• The plane wave expansion for the auxiliary wave functions must be com- 
plete. The plane wave expansion is controlled easily by increasing the 
plane wave cutoff defined as Epw — ^h^G^^^^. Typically we use a plane 
wave cutoff of 30 Ry. 

• The partial wave expansions must be converged. Typically we use one 
or two partial waves per angular momentum (£, to) and site. It should 
be noted that the partial wave expansion is not variational, because the 
partial wave expansion changes the total energy functional and not only 
the basis set. 

We do not discuss here numerical approximations such as the choice of the radial 
grid, since those are easily controlled. 

We mentioned earlier that the pseudopotential approach can be derived as 
a well defined approximation from the PAW method: The augmentation part 
AE = — is a functional of the one-center density matrix Dij defined in 
Eq. E3. The pseudopotential approach can be recovered if we truncate a Taylor 
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expansion of AE about the atomic density matrix after the hnear term. The 
term linear to Di^ is the energy related to the nonlocal pseudopotential. 

r)/\ K 

= E,,if + J2 fn{^n\Vra\^n) + 0{D,^, ~ D^^^f (32) 

n 

Thus we can look at the PAW method also as a pseudopotential method with 
a pseudopotential that adapts to the instantaneous electronic environment, be- 
cause the explicit nonlinear dependence of the total energy on the one-center 
density matrix is properly taken into account. 

What are the main advantages of the PAW method compared to the pseu- 
dopotential approach? 

Firstly all errors can be systematically controlled so that there are no trans- 



ferability errors. As shown by Watson | Watson ct al 1998| andKresse [Kresse et al 1999 



most pseudopotentials fail for high spin atoms such as Cr. While it is probably 
true that pseudopotentials can be constructed that cope even with this situa- 
tion, a failure can not be known beforehand, so that some empiricism remains 
in practice: A pseudopotential constructed from an isolated atom is not guar- 
anteed to be accurate for a molecule. In contrast, the converged results of the 
PAW method do not depend on a reference system such as an isolated atom, 
because it uses the full density and potential. 

The PAW method provides access to the full charge and spin density, which 
is relevant for hyperfine parameters. Hyperfine parameters are sensitive probes 
of the electron density near the nucleus. In many situations they are the only 
information available that allows to deduce atomic structure and chemical envi- 
ronment of an atom. There are reconstruction techniques for the pseudopoten- 



tial approach, which however, are poor man's versions [ Van de Walle et al 1993 1 
of the PAW method. 

The plane wave convergence is more rapid than in norm-conserving pseu- 
dopotentials and should in principle be equivalent to that of ultra-soft pseudo- 



potentials |Vandcrbilt 1990 1. Compared to the ultra-soft pseudopotentials, how- 
ever, the PAW method has the advantage that the total energy expression is 
less complex and therefore is expected to be more efficient. 

The construction of pseudopotentials requires to determine a number of pa- 
rameters. As they influence the results, their choice is critical. Also the PAW 
methods provides some flexibility in the choice of auxiliary partial waves. How- 
ever, this choice does not influence the converged results. 



7 Recent Developments 

Since the first implementation of the PAW method in the CP-PAW code, a num- 
ber of groups have adopted the PAW method. The second implementation was 



done by the group of Holzwarth [ Holzwarth et al 1997 1. The resulting PWPAW 



14 



code is freely available jTackett at al 2001 1. This code is also used as a basis for 



the PAW implementation in the Ablnit project | AbInit|. An inde pendent PAW 
code has been developed by Valiev and Weare [ Valiev et al 1999|. Recently th e 
PAW method has been implemented into the VASP code |Kresse et al 1999 1. 
The PAW method has also been implemented by W. Kromen into the EST- 
CoMPP code of Bliigel and Schroder pliigel et al 2001 1. 

Another branch of methods uses the reconstruction of the PAW method, 
without taking into account the full wave functions in the self-consistency. Fol- 
lowing chemist notation this approach could be termed "post-pseudopotential 
PAW". This development began with the evaluation for hyperfine parame- 
ters from a pseudopotent ial calculation using the PAW reconstruction operator 
[ Van de Walle et al 1993 and is now used in the pseudopotenital approach to 
calculate properties that require the correct wave functions 

The implementation by Kresse and Joubert [ Kressc et al 1999 has been par- 
ticularly useful as they had an implementation of PAW in the same code as the 
ultra-soft pseudopotentials, so that they could critically compare the two ap- 
proaches with each other and LAPW calculations. Their conclusion is that both 
methods compare well in most cases, but they found that magnetic energies are 
seriously - by a factor two - in error in the pseudopotential approach, while 
the results of the PAW method were in line with other all-electron calculations 
using the linear augmented plane wave method. As a short note, Kresse and 
Joubert incorrectly claim that their implementation is superior as it includes 
a term that is analogous to the non-linear core correction of pseudopotentials 
[ Louie ct al 1982| : this term however is already included in the original version 
in the form of the pseudized core density. 

Several extensions of the PAW have been done in the recent years: For ap- 
plications in chemistry truly isolated systems are often of great interest. As 
any plane-wave based method introduces periodic images, the electrostatic in- 
teraction between these images can cause serious errors. The problem has been 
solved by mapping the charge density onto a point charge model, so that the 
electrostatic interaction could be subtracted out in a self-consistent manner 
[ Blochl 1995 [. In order to include the influence of the environment, the latter 
was simulated by simpler force fields using the molecular-mechanics-quantum- 
mechanics (QM-MM) approach |Woo ct al 200C [. 

In order to overcome the limitations of the density functional theory several 
extensions have been performed. Bcngonc [ Bengone et al 2000 implemented 
the LDA+U approach [Anisimov et al 1991 into the CP-PAW code. Soon after 
this, Arnaud [Arnaud ct al 2000 accomplished the implementation of the GW 



approximation into the CP-PAW code. The VASP-version of PAW [ Hobbs et al 2000 



and the CP-PAW code have now been extended to include a noncoUinear descrip- 
tion of the magnetic moments. In a non-coUincar description the Schrodinger 
equation is replaced by the Pauli equation with two-component spinor wave 
functions 

The PAW method has proven useful to evaluate electric field gradients 



[Pctrilli et al 1998 and magnetic hyperfine parameters with high accuracy | Blochl ct al 2000 
Invaluable will be the prediction of NMR chemical shifts using the GIPAW 
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method of Pickard and Mauri jPickard ct al 2001], which is based on their earher 



work [Mauri ct al 1996|. While the GIPAW is implemented in a post-pseudo- 
potential manner, the extension to a self-consistent PAW calculation should be 
straightforward. An post-pseudopotential approach has also been used to evalu- 



ate core level spectra [ Jayawardane et al 2001 1 and momentum matrix elements 



[Kageshima et al 1997 1. 
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